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Abstract 

The IRAS survey data can be used successfully to produce images of extended objects. The 
major difficulties, viz. non-uniform sampling, different response functions for each detector, and 
varying signal- to- noise levels for each detector for each scan, have been resolved. The results 
of three different image construction techniques are compared: co-addition, constrained least- 
squares, and maximum entropy. The maximum entropy result is superior. We present an image 
of the galaxy M51 with an average spatial resolution of 45 arc seconds, using 60 micron survey 
data. This exceeds the telescope diffraction limit of 1 minute of arc, at this wavelength. Data 
fusion is a proposed method for combining data from different instruments, with different spatial 
resolutions, at different wavelengths. Direct estimates of the physical parameters, temperature, 
density and composition, can be made from the data without prior image (re-)construction. An 
increase in the accuracy of these parameters is expected as the result of this more systematic 
approach. 


1 Introduction 

The Infrared Astronomical Satellite (IRAS) surveyed about 95% of the sky, in four broad spectral 
bands centred on 12, 25, 60, and 100 microns, during a ten month period in 1983. Precession at a 
rate of about 1° per day, kept the orbit of the spacecraft remaining perpendicular to the earth-sun 
vector (Figure 1). A semi-overlapping scan strategy was used for the ‘all-sky’ survey. Redundant 
coverage on the time scale of hours was provided by advancing the instrument in elongation by 
half of the width of the focal plane on a subsequent scan, usually the next orbit. For each spectral 
band there were two detector arrays (Figure 2). The arrays were arranged such that the second one 
scanned the same area of sky some 5 to 10 seconds later than the first. Different scans over the same 
area usually intersect at an angles due to precession of the satellite orbit. In addition, most of the 
scans were taken along small circles. Therefore, even small areas of the sky can be very unevenly 
covered . 

•Talk presented at the Workshop on Restoration of IIST Images and Spectra, August 21-22 1990, Baltimore. 
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The nature of the IRAS data is a collection of detector scans cross-cutting the sky at various 
angles. Consequently, the data are not in the form of an image and require non-traditional methods 
of image reconstruction. The term ‘image reconstruction ’ actually does not apply; the procedure 
should be called image construction, because the true scene is unknown. Traditional reconstruction 
routines start with images on an evenly spaced rectilinear grid, of objects convolved with a single 
point-spread function. 

Each of the 59 active IRAS detectors have different Response Functions (RFs). Most detectors 
have rectangular apertures of 0.75-3 (in the in-scan direction) to 4.5-5 arc minutes (in the cross-scan 
direction). Consequently, the spatial resolution is different in the two directions. It is possible to 
improve on the spatial resolution in the cross-scan direction because of the confirmation strategy, 
and because the two rows of detectors for each wavelength band are shifted by half a detector length 
in the cross-scan direction. 

The detector outputs of two consecutive scans over the galaxy M51 at 60 microns wavelength 
are shown in Figure 3. This figure can be regarded as ruled-surface plots of the area. The signal 
is shifted between the two plots due to the half focal plane offset between scans. The main galaxy 
and its satellite NGC 5195 are clearly resolved. Figure 4 shows the centres of the individual sample 
positions of the M51 area, together with the outline of a standard 60 microns detector. 

2 Data representation 

Suppose datum d n is the calibrated value of the n-th sample, taken by detector number i with its 
centre at the position ( x n ,y n ). Detector i has response function Ri(x,y). Note that in general the 
RF profile is rotated according to the scan angle, as in Figure 4. Assume that for each sample this 
rotation has been taken into account in a temporary re-definition of Ri(x,y ), now with axes parallel 
to the axes of the desired map. The measured datum is now the result of a two dimensional integral 
of the sky brightness b(x,y) and Ri(x,y): 

d n = J dxdyRi(x n - x,y n - y)b(x,y) + n„, (1) 

with n n the noise in this n-th datum. Strictly, this is not a convolution since d n is a single number, 
and not a function in the continuous variables x and y\ therefore it is called a sample. 

Digitizing the brightness distribution in pixels reduces the integral to a summation: 

M 

d n — ^ ^ r nm b m T n n , (2) 

m= 1 

with 

r nm = / dxdyRi(x n - x,y„ - y). (3) 

J area pixel m 

Equation 2 can be interpreted as one equation in M unknown discrete brightnesses b mi with known 
coefficients r nm . There are in total N samples falling (partly) inside the map, forming a set of N 
equations in M unknowns (e.g. in Figure 4 N ~ 700). This can be written as the matrix equation 

d — Rnm Hn, (4) 
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with d = (cfi, *■ • , dtf) and n = (ni, • • ■ ,un), the data and the noise vector, respectively. The matrix 
Rnm is the N x M response matrix, in which the nth row is the set of coefficients R nm for datum 
d n . The unknown sky is represented by b = (6j, • • • ,6 m), a vector of length Af, obtained through 
stacking of the rows (or columns) of the desired map. The RFs are normalized to unit volume. 

For samples which overlap the map boundary, it is assumed that the mean intensity just outside 
the map is the same as just inside the boundary. In this way boundary effects can be apodized. 

Formally, the standard deviation n n of each datum d n is separately represented in Equation 4. 
The noise level is estimated from the available data, viz. signal plus noise, using a zero-sum filter. 
Differences in noise level of a factor three for the same detector, in consecutive orbits have been 
measured. 

Summarizing, the non-uniform sampling, the different RF for each detector and their arbitrary 
rotations are incorporated in the response matrix. The different noise levels per detector and per 
scan are represented in the noise vector. The solution b is defined on a regular pixel grid. Since the 
individual response functions are used, no special provisions are necessary for the small detectors. 
The image construction has become a numerical mathematical problem which can be solved by many 
different methods. 

Three methods have been compared, co-addition, constrained least-squares, and maximum en- 
tropy (see Bontekoe et al., 1991). The best results are obtained with maximum entropy, and are 
presented here. 


3 Image construction 

Even when the response matrix Rnm is perfectly known, the recovery of the original scene is 
mathematically impossible. Given the data d, the solution b is not unique. The problem is called ill- 
posed. In addition, the image (re-)construction problem is also ill-conditioned. The ‘signal-to-noise’ 
ratio of the final map is usually orders of magnitude worse than of the input data. Nevertheless, 
astronomically meaningful images can be obtained through regularization of the problem; maximum 
entropy is such a regularization. 

In addition to solving for b the standard deviation is also computed for each pixel, and serves 
as an error map. Although this ignores the covariance of neighbouring pixels, we use this error map 
to give us some indication of the reliability of features in the solution b . 

The images are constructed on a grid of 60x60 pixels, for the same map area as Figure 4 (15 x 15 
arc minutes). Since there are N ~ 700 samples and M = 3600 unknowns, the system of equations 
is underdetermined. 

The MEMSYS3 software package of Skilling (1989) and Gull (1989) is used to obtain a maximum 
entropy estimate b from M51 data. Until recently, the common procedure was to maximize the 
entropy under the constraint that the goodness- of- fit statistic equals the number of data points. This 
approach, however, has two major drawbacks. First, within this frame-work there is no consistent 
method to estimate the error map cr^, although the maximum entropy solution b can be found 
straightforwardly. Second, the \ 2 = N criterion gives no allowance for the fact that significant data 
will generate structure in the image fc. This structure is like a set of parameters being fitted from 
the data. The effective number of parameters G supporting the underlying the structure should be 
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subtracted from N to obtain the number of degrees of freedom with which to compare x 2 (see also 
the discussion in Gull, 1989). 

The Bayesian estimate for the present case is G ~ 400. The remaining N — G ~ 300 is the 
number of degrees of freedom for this problem, and this should be used for the \ 2 statistic. 

The maximum entropy solution shows fine details, especially in areas of low brightness where 
the other methods fail. Although such maps can be used sensibly, inferences ought to be made by 
overlaying them with masks and computing integrals over the masked solution. 

The Figures 5 and 6 are the result of coarsening the maximum entropy solution, originally 
computed on 60x60 pixels, to an effective 30x30 grid by using 900 2x2 pixels boxcar masks, being 
1 in the square of interest and 0 outside. The 60x60 grid is retained, however, and each group of 
four pixels is assigned the same value b. Finally, this coarsened map is smoothed again using a 2x2 
boxcar filter. 

Division of brightness map by the error map yields Figure 6. 

Most of the spiral structure in Figure 5 coincides well with the spiral arms of the Ha image at 8 
arc second resolution in van der Hulst et al. (1988). Extensions of the spiral arms at low brightness 
levels in the maximum entropy solution line up very well with outer parts of the spiral arms in Ha. 
The brightest point does not coincide with the nucleus in the Ha, but lies about 25 arc seconds 
to the South. The second source in the centre also has no visual counterpart and coincides with 
an inter-arm region. The companion galaxy seems resolved into a strong point source, towards its 
South-East, which coincides exactly with a sharp maximum in Ha and extended structure towards 
the North and West. 

Figure 6 summarizes the maximum entropy result in a statistical sense. The peaks in the galaxy 
and companion are 8 cr and 9(7 detections, respectively. The spiral arms are 1-3 a detections, and 
the secondary peak in the nucleus of M51 a 3cr detection. A unidentified source near the Northern 
boundary of the map is a 4<7 detection, although its brightness is low. This might be an artifact 
from the treatment of the boundary. 

A difficult issue is the final spatial resolution in the map. Since MEMSYS3 finds G ~ 400, it 
is tempting to distribute this number evenly over the image as the number of independent picture 
elements, yielding an average spatial resolution of 45 seconds of arc. This is an improvement above 
the diffraction limit of the telescope as well as the classical limit imposed by sampling theory, both 
of which are 1 arc minute. Spatial resolution, however, is dependent on the ‘signal-to-noise’ in the 
original data and consequently non-uniform over the map. Therefore the 400 fitted parameters can 
not be evenly distributed, and areas of high brightness can have a better spatial resolution than 
average. The reverse is true for low brightness areas. 

Although the IRAS survey mission was not intended to produce images, the major difficulties, viz. 
non-uniform sampling, different RFs for each detector, and different signal-to-noise levels for each 
detector for each scan, have been dealt with. Astronomically meaningful images can be produced 
of the far-infrared sky from the IRAS survey data, but they have to be constructed by advanced 
numerical techniques. Overall, the MEMSYS3 result is superior to the others. The images produced 
by MEMSYS3 show better spatial resolution, are non-negative and show plausible structure even at 
low brightness levels. 
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4 Data Fusion 


One of the centra] questions in astrophysics is the correct identification of radiative processes op- 
erating in a given source. A fresh attack is proposed on the classification of physical processes in 
a complex field from a fundamentally new direction. The goal is to take data of a given field, ob- 
tained at a variety of wavelengths and spatial resolutions, and produce images of operative physical 
processes and the corresponding parameters. 

First the data acquisition process is briefly described. Emphasis is given to the correct dimensions 
of all quantities. The sky brightness B(8, <f>, A) (in Wm~ 2 m- l sr~ l ) is a function of position of the sky 
and of wavelength A. When the instrument is pointed in the direction (0;,<£,), the response 
function, representing the blur, is defined R(6, <j> t 6i y <f>i) (dimensionless). In addition the blurred 
signal must pass a colour filter F( A) (dimensionless). Each datum of measured flux dj (in Wm~ 2 ) 
is now the result of both ‘convolutions’, 


d, i = J J d0d<t>R(0,4>,6i,<t>i) j d\F(\)B(0,<j>,\). 


(5) 


The important fact to notice is that the data acquisition is assumed to be linear with the input 
signal, which is the essential task for the calibration. It is this linearity that allows decomposition 
of the sky brightness into various components. Note also the difference between B here and b of the 
previous sections; the latter is the brightness integrated over the colour filter. The effect of noise is 
well understood in linear problems, and does not affect the theoretical analysis. 

The ensemble of number densities n 8 (r y 0 y <j> y T s ) (in m~ 3 ) of sources, indexed with s y can be 
estimated, when the source’s radii R s (in m), and template spectra / a (A,T s ) (in Wm~ 2 m~ 1 $r~ 1 ) 
are given. The best fitting temperature T s (in degrees K) of the various n 8 results as part of the 
solution. An infinitesimal volume dv = r 2 drd£l (in m 3 ), radiates with a spectral power per unit 
wavelength dw (in Pfm" 1 ) equal to 


dw(r,0J) = dv n,(r,0,<j>,T,)A^R]l,(\,T,), 
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( 6 ) 


assuming isotropic radiation by all sources. The contribution to the sky brightness of dv is 


dB(0J, A) 


MbM) -r(r, A) 
47 rr 2 dfi 


(7) 


where r(r, A) is the optical depth (dimensionless). Integrating over the line of sight, the sky bright- 
ness becomes 


8(0,4, A) = / dre~ r ^-^ V n s (r,0,<j>,T,) R 2 , I S (\,T 3 ). 
' J 0 . 


( 8 ) 


Under the condition that the absorption, represented by r, and the emission I s are independent of 
the local radiation field, the sky brightness B is a linear function of the densities n s . 

Since the data acquisition is assumed linear with brightness, there is linear relation between 
the data d{ and n 8 the physical parameters of the distribution of matter in space. It requires 
simultaneous solution of the equation of radiative transfer and the instrumental inversion. However 
complicated, the problem is linear, and linear problems can be solved by many numerical techniques. 
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The critical factor in resolving three-dimensional structures in the universe is knowledge of the factor 
T (r,0 } <f>, A). The full solution n,(r,0,<£, T,) for all sources describes the composition, temperature, 
and density of all visible matter. If this information is not present in the data, upper bounds for the 
desired densities are still a useful result. 

If one is less confident in the quality of the data or ones knowledge of the absorption, one can 
apply a two dimensional version of the theory in which only projected densities are defined. 
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Figure captions 

Figure 1: Schematic drawing of IRAS orbital geometry. 

Figure 2: IRAS focal plane. Of the 62 infrared detectors the 3 filled-in were inoperative; the 
cross-hatched detectors showed a higher noise level. The normal scan direction of images is shown. 

Figure 3: Detector scans from two consecutive scans covering the galaxy M51 (60 microns). The 
scan direction is from right to left and each has a length of 0.5 degree. The detector scans are 
displaced vertically corresponding to their cross-scan position in the focal plane. The figures can be 
regarded as ruled surface plots. The maximum signal in both plots is 430 detector units. 

Figure J: Positions of the samples in the area of M51 (60 microns). The circles represent the 
two smallest detectors in the band, viz. detectors 11 and 31 (see Fig. 1). The scans roughly run 
from top-left to bottom-right. The area is covered by about 700 samples. The fat contours outline 
the RF profile of a normal size detector, at the 90, 50, 10, 2, and 1% level of its maximum. The 
detector is centred at the position indicated with an asterisk. A grid with a one arc minute spacing 
is superimposed. 

Figure 5: Map of M51 (60 microns) from the MEMSYS3 maximum entropy method. Area is the 
same as in Figure 4. The lowest contour is at a level 150; subsequent contours with increments of a 
factor 2. 

Figure 6. Signal- to- noise map of the MEMSYS3 solution of M51. Lowest contour is at 1 ( 7 , 
subsequent contours are separated by 1 a. 
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